Timescale mediates the effects of environmental controls on water temperature in mid- to low-order streams

Adequate management and conservation of instream thermal habitats requires an understanding of the control that different landscape features exert on water temperatures. Previous studies have extensively explored the influence of spatial scale on these relationships. However, the effect of temporal scale remains poorly understood. Here, we use paired air–water mean daily and monthly summer temperatures collected over four years from 130 monitoring stations in Japanese mid- to low-order streams to investigate whether perceived effects of different environmental controls on water temperature are dependent on the timescale of the temperature data, and whether those dependencies are related to the spatial scale at which these controls operate. We found a clear pattern for the significant cooling effect, high relative importance and strong dominance exerted by the riparian forest cover on daily temperatures at the reach scale becoming dampened by concomitant increases associated to the proportion of volcanic geology on monthly temperatures at the catchment scale. These results highlight the importance of contextualizing the effects of environmental controls on water temperatures to the timescale of the analysis. Such dependencies are particularly important for the management and conservation of instream thermal habitats in a rapidly warming world.

Increasing thermal stress resulting from human disturbances poses a risk to freshwater biota, given the control of temperature on physiological processes 1 , prompting them to respond in situ by exploiting local thermal refugia 2 or by relocating to colder tributaries and upstream reaches 3 . As a result, characterizing instream thermal habitats and their controlling factors at both local and broader, basin-wide scales is becoming increasingly important for the integrated management and conservation of fluvial ecosystems and associated biota 4,5 .
Multiple climatological, hydrological, morphological and geological factors influence stream water temperature across spatial scales through their effect on heat exchange processes occurring at the stream surface and streambed interfaces 6 . A growing body of literature demonstrates how spatial heterogeneity in the characteristics of river networks and their catchments can shape localized responses of stream temperature sensitivity, which can vary significantly across space even in regions that are geographically close and climatically similar [7][8][9] . This makes the extrapolation of simple air-water temperature relationships beyond monitored sites across river networks problematic. Yet predictions of water temperature are needed at scales meaningful to management (i.e., catchment scale). Recent modelling efforts try to circumvent this issue by explicitly incorporating the effect of landscape and geomorphological characteristics into the prediction of stream water temperatures [10][11][12] .
Thermal regimes are also geographically distinct and scale dependent in the time domain 13 . Previous studies have shown that air-water temperature relationships vary significantly with the temporal resolution of the data, where higher slopes and lower intercepts are usually associated when regressing water on air temperature at increasing timescales 6 . However, much less is known about how temporal scale can modify perceived effects of other environmental controls on stream water temperatures over and above the effect of air temperatures. Using spatial stream network models, Steel et al. 14 found that relationships between landscape predictors and water temperature in the Snoqualmie River changed among temperature metrics defined for different temporal windows and timescales. For example, while increasing elevation correlated strongly with decreasing summer mean temperatures, this relationship was reversed for mean daily thermal range and lost for temperature variability   S2). This study is based on data collected during the summer (June-September) between 2018 and 2021 from a total of 130 paired (air-water) temperature monitoring stations (Fig. 1). Stations were allocated across each river basin covering as much of their environmental gradients as possible subject to logistic and accessibility constraints. The monitored network covers only low-to mid-order streams (Strahler order 1 to 4) due to operational constraints and the high level of human alteration existing in the larger, non-wadable rivers of these four catchments.
Air and water temperatures were recorded at each site using Onset Hobo UA-001-64 (accuracy ± 0.53 °C) and Gemini TG-4100 (± 0.5 °C) temperature data loggers. Water loggers were installed in non-turbulent flowing sections of the stream housed in short sections of perforated PVC pipes to facilitate water exchange and shield them from direct sunlight. Air loggers on perforated PVC pipes were attached to nearby trees (~ 2 m high) in shaded areas by the stream bank 24,25 . Loggers were set to record temperatures every hour, which were then converted into mean daily and monthly temperatures using the 'timeAverage' function of the R package ' openair' 26 at a 0.75 unit threshold. In other words, a given day or month was computed from the corresponding hourly data if at least 75% of the hours for that temporal unit were available. Otherwise the unit was taken as a missing observation. Visual inspection of the temperature series suggested some suspected periods of air exposure of water loggers in a few sites (Fig. S4), which were subsequently treated as missing observations. Missing observations occurred also sporadically due to issues such as faulty or lost loggers. Nonetheless, the proportion of missing observations within individual series was small across sites (1.79 ± 5.79%), resulting in a final total of 38,798 daily and 1324 monthly valid paired temperature observations with an average effective series duration by site of 10.19 ± 2.46 months (234.85 ± 144.23 days).
Landscape variables. River networks and catchments were derived from the conditioned 10-m grid digital elevation model (DEM) form the Geospatial Information Authority of Japan (available from https:// fgd. gsi. go. jp/ downl oad/) based on flow accumulating and minimum length thresholds of 5000 (~ 0.5 km 2 minimum contributing catchment) and 100 cells (~ 1 km), respectively, using the Spatial Analyst Tools in ArcGIS 10.7.1. River network topology and different environmental layers were then queried at each monitored site to extract our set of predictor variables (Table 1).
Statistical models. We used linear mixed effects models (LMMs) for prediction of mean daily and monthly water temperatures following the general notation: where Tw ij and Ta ij are the i-th observation of mean water and air temperatures for site j; α is the (grand) intercept; βX i,j represents the corresponding linear combination of site-specific environmental covariates and their corresponding regression coefficients; ϵ i,j are the observation-specific residuals; υ j,k is a random effect that includes a random slope and intercept for air temperature at the site level, accounting for the different air-water temperature relationships observed across sites (Fig. S5); and h(s j , ϕ i ) is a continuous temporal autoregressive (1) Tw i,j = α + β air Ta i,j + βX i,j + ϑ j,k + h s j , ϕ i + ǫ i,j Table 1. Description of catchment-and reach-scale environmental covariates measured at each monitored site including related hydrological and stream heat exchange processes (+ /− symbols indicate expected positive/ negative relationships), and the source and format of the data sets used in this study to derive each variable. www.nature.com/scientificreports/ (CAR) structure of order 1 at site level used in the daily model (see below) to account for temporal correlation in model residuals defined by φ, the correlation between two observations one unit of time apart, and s as the time covariate grouped by site 27 . We initially considered the inclusion of catchment as the highest level of a hierarchical nested random effect structure. However, we decided to include it as a fixed-effect categorical variable due to the limited sample size at the higher hierarchical level (i.e., 4 catchments). A bare minimum of five groups is usually considered necessary at the highest hierarchical random effect level to make meaningful estimates of the variation among random effects (group-level) in hierarchical mixed-effects models 28,29 . Examination of the normalized residuals by site for the daily model showed a clear, consistent pattern of a gradually decaying autocorrelation function and a partial autocorrelation function with a strong significant correlation (~ 0.8-0.9) at lag 1 dropping sharply below the significance threshold thereafter (Fig. S6). Therefore, we included a continuous autoregressive structure of lag 1 into the daily models, which flexibly allows for the existence of gaps in the time series of records 27 . These temporal autocorrelation patterns were absent in the monthly model (Fig. S6). Therefore, the temporal correlation structure was not added to the monthly model.
Despite the limited number of observations at the site level for the monthly model (10.19 ± 2.46 observations), simulation studies demonstrate that, when the number of groups is sufficiently large (i.e., 130 sites in our case), neither fixed nor random effects estimates are affected by small sample size within groups even in extreme situations 30,31 . Sample size at site level for the monthly model was hence considered to be sufficient to obtain unbiased estimates of the fixed effects for the monthly model at the population level (i.e., across sites), which constitute the focus of this study.
Model performance was assessed by mean of the conditional and marginal coefficients of determination R 2[,32,33 , representing respectively the total variance explained by the model (fixed and random components together) and the variance explained by the fixed effects. The root mean square error (RMSE) was also computed for each model as a measure of goodness-of-fit. All metrics were computed using the function 'model_performance' from the R package 'performance' 34 . Covariates were assessed for multicollinearity on the fitted models using variance inflation factors (VIFs) with a cut-off value of 10 35 . All predictors had VIFs well below that threshold in both models.

Relative importance of predictors and dominance analysis.
We assessed the relative importance of each predictor in the models using Dominance Analysis (DA) 19,36 . DA is a technique to rank-order the predictors in a model in terms of relative importance by comparing the additional contributions from each predictor to the variance explained by all subset models containing all possible combinations of the rest of the predictors. The additional contribution of a predictor to a given subset model is measured as the increase in the variance accounted for by that model after the addition of that predictor. Dominance can then be established among any pair of predictors present in the original model, where one predictor is said to dominate over another if its additional contribution is greater than that of the other predictor over all subset models not including any of the two predictors 19 . Azen and Budescu 36 established three types of dominance that are hierarchically related. At the highest level, complete dominance of a predictor over another is established if the additional contribution of that predictor is consistently greater than the contribution of the other predictor over all possible subset models that do not include any of the two predictors. When this is not the case, weaker levels of dominance can still be established between pairs of predictors. Conditional dominance occurs when the average additional contribution of a predictor for all subset models within each model size is greater than that of another predictor. Finally, general dominance, representing the weakest level in the dominance hierarchy, occurs when the additional contribution of a predictor averaged across all model sizes is greater than that of the other predictor. Although several other related methods for assessing the relative importance of predictors based on an 'all subsets' approach exist 37 , their estimations focus only on average contributions over all predictor orderings (i.e., general dominance in DA). Therefore, DA provides the additional important insight given by the other relative importance measures (i.e., complete and conditional dominance) 38 .
Given our interest is in assessing the control of environmental predictors on water temperatures, we used a model fitting water temperatures to air temperatures plus the catchment effect as a baseline against which relative contributions of the environmental predictors were evaluated in terms of changes in marginal R 2 while keeping the random component of the models invariant. This allowed us to assess predictor contribution to explain variability in water temperatures after accounting for the effect of air temperatures. All statistical computations were done in R version 3.6.3 39 .

Results and discussion
The proportion of variance in water temperatures explained by the models increased with coarser temporal resolution with marginal and conditional R 2 of 0.71/0.81 and 0.77/0.93 for mean daily and monthly temperatures, respectively. RMSEs (population-level) showed a reversed trend, with mean deviations of predicted water temperatures from observed water temperatures decreasing from 1.83 °C for mean daily temperatures to 1.41 °C for mean monthly temperatures. Models at both timescales improved the performance of the baseline model (R 2 day = 0.5, RMSE day = 2.39 °C, R 2 month = 0.57, RMSE month = 1.94 °C) significantly (p < 0.0001 ANOVA based on Maximum Likelihood model estimates), demonstrating the utility of the inclusion of environmental variables to capture variability in air-water temperature sensitivity across sites.
Air temperature had a highly significant positive effect on water temperature with a regression coefficient twice as big in the monthly compared to the daily model (Table 2). Steeper regression slopes and lower intercepts often result when increasing the timescale in simple water-air linear regression models because averaging over longer timescales decreases the variance in the data by reducing the effect of the time lag or transient response of stream water temperatures to changing air temperatures, typically operating on the order of hours to days 40 www.nature.com/scientificreports/ catchment effect was also significant in both models (Table 2); the southern catchments (Kiso and Hiji) having higher water temperatures on average relative to the northern (Sorachi and Teshio) catchments. Some environmental predictors had consistent significant effects on water temperature irrespective of the timescale of the temperature data ( Table 2). Increasing proportion of quaternary volcanic rocks within the contributing catchment and higher reach elevation had a highly significant negative effect on water temperatures while increasing proportion of cultivated land in the catchment had a significant positive effect. The former two variables had also the highest importance in terms of relative contribution to the variance explained by the fixedeffect component of the model relative to the baseline (Fig. 2a-b), with the proportion of quaternary volcanic rocks accounting for a 31.48% and 38.41%, and reach elevation a 27.06% and 26.53% of the total variance in the daily and monthly models.
Other predictors, however, showed clear differences in their effects and relative importance between the two models ( Table 2, Fig. 2a-b). At the reach scale, the cooling effect of riparian cover became weaker with increasing timescale, its effect size halved and shifting from highly significant to marginally non-significant in the monthly model (Table 2); its relative importance reduced from 23.52 to 16.93% (Fig. 2a-b). On the other hand, at the catchment scale, increasing proportion of pre-quaternary extrusive volcanic formations had a significant negative effect on monthly water temperatures that was absent from the daily model (Table 2), while the relative importance of pre-quaternary and quaternary volcanic rocks increased by 4.65%, and 6.9% in the monthly relative to the daily model ( Fig. 2a-b).
These results were further supported by dominance analysis (Fig. 2c-d). Despite being the most dominant predictor in both models, the proportion of quaternary volcanic rocks experienced a strong weakening of dominance links on other predictors between timescales, changing from exerting complete dominance over all other predictors but reach elevation in the monthly model to having only general dominance (i.e., lowest hierarchical level of dominance) over riparian forest cover and conditional dominance over the proportion of catchment cultivated land in the daily model. On the other hand, riparian forest cover and cultivated land use increased both one level in dominance hierarchy over the proportion of pre-quaternary volcanic rocks in the daily relative to the monthly model (Fig. 2c-d).
Together, notwithstanding the fact that reach-and catchment-scale covariates remained important at both temporal resolutions, these results point towards a clear temporal dependency of the relative importance and dominance of environmental controls on stream water temperature. The stronger effects, higher relative importance and dominance exerted by reach-scale predictors related to controls governing incident solar radiation and local heat exchange processes at the stream surface observed at the daily timescale became dampened at the monthly timescale by predictors related to catchment-scale hydrological and hydrogeological processes controlling the partitioning, storage and routing of water into the streams. These results hint at potential interdependencies between temporal and spatial scales that will be explored in more detail in future studies. www.nature.com/scientificreports/ Cooler daily summer stream temperatures are often associated to the effect of riparian canopy limiting incoming radiation at the water surface resulting in a significant reduction of the net surface heat flux compared to open or non-vegetated riparian zones [43][44][45] . Similar effects have also been reported, albeit less frequently, in studies focusing on coarser timescales (i.e., monthly or annual temperatures) 46 , perhaps because at these temporal resolutions the interest is often on larger spatial scales where the focus is on predictors related to catchment-scale processes such as land use cover 47 . In any case, our results suggest that riparian cover exerts a larger relative control on water temperature at finer than coarser timescales. Riparian vegetation can influence short-term (e.g., daily) variation in water temperatures via lateral flow of surface and ground water into the stream channel through interception, infiltration and soil saturation processes, particularly in small, steep headwater streams 48,49 .
The relative importance and dominance niche left by riparian cover in the monthly model was filled by the proportion of solid volcanic geology in the contributing catchment. Although the location, direction and rate of groundwater-surface water interactions are controlled by multiple processes that operate over different spatial and temporal scales 50 , basic information on solid geology at catchment level can be an effective proxy for groundwater-streambed interactions 51 . Together with the type of formation, geological age is important as it relates to the processes of weathering, soil development and formation of secondary porosity in underlying rocks, all of which are associated to aquifer recharge and surface-groundwater dynamics 20 . Young extrusive volcanic formations often have high porosity and permeability offering deep, stable aquifers. In our study catchments, young volcanic formations comprised primarily mafic basaltic and andesite rocks, both formations capable of forming prolific unconfined aquifers in which groundwater flow freely through vesicular and fracture systems 20 . On the other hand, older extrusive volcanic rocks often experience a decreasing trend in porosity and hydraulic conductivity with increasing geological age resulting from the filling of vesicles and fractures by secondary material. In such settings, groundwater is largely stored in less stable, shallow aquifers with short flow paths to streams. A study conducted at 70 gauging stations across Japan, including sites located in three of our four study catchments, found a strong correlation between the type and age of dominant catchment geology and base flows with catchments dominated by quaternary volcanic formations supporting the most abundant flows whereas those under older volcanic formations yielded intermediate flows 21 . Shorter residence times and shallower storage also mean www.nature.com/scientificreports/ warmer groundwater temperature and greater sensitivity to seasonal and long-term changes in air temperatures 52 . The observed weaker but significant negative effect of old volcanic rocks on monthly water temperatures likely reflects to some extent these contrasting conditions. Catchment cultivated land use had a consistent and significant warming effect on water temperatures irrespective of resolution and, together with riparian forest cover, increased in relative importance and dominance in the daily model. Agricultural land use is often associated to increased evapotranspiration and surface runoff and decreased base flow and water yield 53 ; all factors that can contribute towards warming of thermal regimes in rivers particularly during dry seasons 54,55 . Decreased time lag responses between precipitation and streamflow associated to agricultural use 56 may have contributed towards the observed higher signature of this predictor on mean daily relative to monthly water temperatures.

Conclusions
Our perception of nature, intrinsically dynamic in time and space, is inevitably conditioned by the scales and degree of detail we observed it with. Paraphrasing Simon A. Levin "there is no single natural scale at which ecological phenomena should be studied; systems generally show characteristic variability on a range of spatial, temporal, and organizational scales" 57 . To understand ecological patterns, and the processes that produce them, it is necessary to analyze and compare them across multiple spatial and temporal scales. The role of scale in defining patterns and processes is especially important to landscape ecology because landscapes comprise many heterogeneous components that interact and show dependence across spatial and temporal scales 58,59 .
Previous studies have documented complex temporal and spatial dependencies in the contribution of climatological, topographical and geological controls to hydrological and hydrochemical responses of streams 14,60,61 . For example, Karlsen et al 62 found that catchment soil characteristics correlated significantly with specific discharge in a boreal catchment at shorter timescales (monthly to daily), but not when discharge was aggregated into annual timescales. In our study, although reach-and catchment-scale covariates remained important in both models, we found a clear signal for the strong control exerted by the riparian forest cover on daily water temperatures to become dampened by the effect of predictors related to catchment geology in the monthly model. Previous studies have documented the strong, direct control exerted by the type and extent of riparian cover on summer stream temperature variability over short (hour to daily) temporal scales through their influence on the amount of incoming solar radiation received at the stream surface 44 . Increasing the time span over which temperatures are averaged will effectively reduce this variance. This could reduce the effect size, importance and dominance of riparian cover in explaining differences in temperatures across sites between the two models relative to other covariates. On the other hand, groundwater-dominated streams, such as those found in catchments dominated by extrusive volcanic formations in our study, are characterized by cooler, stable summer thermal regimes which signature should be little influenced by the effect of averaging temperatures over increasing time scales. This could explain the observed time-dependent trade-offs between the two covariates in terms of their relative contribution towards explained variance in the models.
Climate change is altering the abundance, distribution and phenology of species across land and oceans, reshuffling biodiversity on a global scale, and modifying ecosystems with dire implications to human wellbeing 63,64 . Species may counter the effects of climate change through phenotypic and evolutionary adaptation although evidence suggest such responses are often outpaced and imperfect 65 . Freshwater ecosystems are at the frontline of these changes because the linear, dendritic nature of hydrological networks constrain biota responses and concentrate compounded impacts from climate change and the many human activities associated to the use of water as a resource 66 . Robust prediction of water temperatures at reach scale across river networks is therefore necessary for the development of sound river management practices that are resilient and adaptive to the effects of future climate change. Statistical regression models that incorporate heat exchange processes in rivers through the incorporation of relevant environmental proxies represent a useful tool to achieve this goal in a simple, flexible way 6 . Our results highlight the importance of contextualizing such relationships to the temporal resolution of the temperature data informing the models.

Data availability
All data sets used for generating the environmental predictors and hydrological networks are publicly available as referenced in Table 1. Air and water temperature data generated by this study is available from the corresponding author on reasonable request. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.